//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __math3d_hpp__
#define __math3d_hpp__

#include <math.h>
#include <stdlib.h>

#if defined(_WIN32) || defined(_WINDOWS)
template <class T>
inline long int lround(T x)
{
  if ( x >= 0 )
    return (long int)(x + 0.5);
  else
    return (long int)(x - 0.5);
}
#endif

#if defined(USE_ERRORPREC) || defined(USE_GMP)
  #include <gmp.h>
  #ifndef GMP_MPF_BITS
    #define GMP_MPF_BITS 256
  #endif
#endif

typedef int int32;
typedef unsigned int uint32;
//typedef long long int64;
//typedef unsigned long long uint64;

#if defined(USE_ERRORPREC)
  struct prec_t
  {
    inline prec_t ()
    {
      mpf_init2(v2, GMP_MPF_BITS);
    }

    inline prec_t(double v) : v1(v)
    {
      mpf_init2 (v2, GMP_MPF_BITS);
      mpf_set_d (v2, v);
    }

    inline prec_t(const prec_t& c) : v1(c.v1)
    {
      mpf_init2 (v2, GMP_MPF_BITS);
      mpf_set (v2, c.v2);
    }

    inline prec_t& operator= (const prec_t& c)
    {
      v1 = c.v1;
      mpf_set (v2, c.v2);
      return *this;
    }

    inline prec_t& operator= (double v)
    {
      v1 = v;
      mpf_set_d (v2, v);
      return *this;
    }

    inline ~prec_t ()
    {
      mpf_clear (v2);
    }

    friend inline bool operator== (const prec_t& l, const prec_t& r) { return l.v1 == r.v1; }
    friend inline bool operator!= (const prec_t& l, const prec_t& r) { return l.v1 != r.v1; }
    friend inline bool operator< (const prec_t& l, const prec_t& r) { return l.v1 < r.v1; }
    friend inline bool operator<= (const prec_t& l, const prec_t& r) { return l.v1 <= r.v1; }
    friend inline bool operator> (const prec_t& l, const prec_t& r) { return l.v1 > r.v1; }
    friend inline bool operator>= (const prec_t& l, const prec_t& r) { return l.v1 >= r.v1; }

    inline prec_t& operator += (const prec_t& r) { v1 += r.v1; mpf_add (v2, v2, r.v2); return *this; }
    inline prec_t& operator -= (const prec_t& r) { v1 -= r.v1; mpf_sub (v2, v2, r.v2); return *this; }
    inline prec_t& operator *= (const prec_t& r) { v1 *= r.v1; mpf_mul (v2, v2, r.v2); return *this; }
    inline prec_t& operator /= (const prec_t& r)
    {
      v1 /= r.v1;
      if ( mpf_cmp_d (r.v2, 0) != 0 )
        {
          mpf_div (v2, v2, r.v2);
        }
      else
        {
          mpf_set_d (v2, 0);
        }
      return *this;
    }

    inline prec_t operator- () const
    {
      prec_t ret;
      ret.v1 = -v1;
      mpf_neg (ret.v2, v2);
      return ret;
    }

    inline friend prec_t operator+ (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      ret.v1 = l.v1 + r.v1;
      mpf_add (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator- (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      ret.v1 = l.v1 - r.v1;
      mpf_sub (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator* (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      ret.v1 = l.v1 * r.v1;
      mpf_mul (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator/ (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      ret.v1 = l.v1 / r.v1;
      if ( mpf_cmp_d (r.v2, 0) != 0 )
        {
          mpf_div (ret.v2, l.v2, r.v2);
        }
      else
        {
          mpf_set_d (ret.v2, 0);
        }
      return ret;
    }

    inline prec_t sqrt () const
    {
      prec_t ret;
      ret.v1 = ::sqrt(v1);
      mpf_sqrt (ret.v2, v2);
      return ret;
    }

    inline prec_t abs () const
    {
      prec_t ret;
      ret.v1 = fabs(v1);
      mpf_abs (ret.v2, v2);
      return ret;
    }

    // returns error in percentage
    double error () const
    {
      if ( mpf_cmp_d(v2, 0) == 0 )
          return v1 == 0.0 ? 0 : 9e9;
      mpf_t e;
      mpf_init2 (e, MPF_BITS_USE);
      mpf_set_d (e, v1);
      mpf_div (e, e, v2);
      double ret = fabs(100.0 * mpf_get_d(e) - 100.0);
      mpf_clear(e);
      return ret;
    }

    inline double d() const { return v1; }
    inline double rd() const { return mpf_get_d(v2); }
    inline void correct () { v1 = mpf_get_d(v2); }

  private:
    double v1;
    mpf_t v2;
  };

  inline prec_t prec_t_sqrt (const prec_t& v) { return v.sqrt(); }
  inline double prec2double (const prec_t& p) { return p.d(); }
  inline size_t prec2size_t (const prec_t& p) { return (size_t)p.d(); }

  inline double cos (const prec_t& v) { return cos(v.d()); }
  inline double acos (const prec_t& v) { return acos(v.d()); }
  inline double sin (const prec_t& v) { return sin(v.d()); }
  inline double asin (const prec_t& v) { return asin(v.d()); }
  inline double tan (const prec_t& v) { return tan(v.d()); }
  inline double atan (const prec_t& v) { return atan(v.d()); }
  inline double floor (const prec_t& v) { return floor(v.d()); }
  inline double ceil (const prec_t& v) { return ceil(v.d()); }
  inline double round (const prec_t& v) { return round(v.d()); }
  inline double lround (const prec_t& v) { return lround(v.d()); }
  inline double prec_t_pow (const prec_t& v, const prec_t& f) { return pow(v.d(), f.d()); }
  inline double prec_t_log (const prec_t& v) { return log(v.d()); }
  inline prec_t fabs (const prec_t& v) { return v.abs(); }

  // check for error in value v, print with name, and give an error level to report (percentage)
  inline void do_prec_checkerror(const char *path, int lineno, const prec_t& v, const char *name, prec_t level = 10.0)
  {
    double e = v.error();
    if ( e > level )
        printf ("%s:%u: %s off by %.2e%%\n", path, lineno, name, e); 
  }

  #define prec_checkerror(a,b,c) do_prec_checkerror(__FILE__,__LINE__,a,b,c)
  inline void prec_clearerror(prec_t& v) { v = prec2double(v); }
  inline void prec_correcterror(prec_t& v) { v.correct(); }

  #define PREC_ERROR
#elif defined(USE_FLOAT)
  typedef float prec_t;
  #define prec_t_sqrt sqrtf
  #define prec_t_log logf
  #define prec_t_pow powf
  #define prec2double(v) ((double)(v))
  #define prec2size_t(v) ((size_t)(v))
  #define prec_checkerror(a,b,c)
  #define prec_clearerror(a)
  #define prec_correcterror(a)
  #define PREC_FLOAT
#elif defined(USE_DOUBLEDOUBLE)
  typedef long double prec_t;
  #define prec_t_sqrt sqrtl
  #define prec_t_log logl
  #define prec_t_pow powl
  #define prec2double(v) ((double)(v))
  #define prec2size_t(v) ((size_t)(v))
  #define prec_checkerror(a,b,c)
  #define prec_clearerror(a)
  #define prec_correcterror(a)
  #define PREC_DOUBLEDOUBLE
#elif defined(USE_GMP)
  struct prec_t
  {
    inline prec_t ()
    {
      mpf_init2(v2, GMP_MPF_BITS);
    }

    inline prec_t(double v)
    {
      mpf_init2 (v2, GMP_MPF_BITS);
      mpf_set_d (v2, v);
    }

    inline prec_t(const prec_t& c)
    {
      mpf_init2 (v2, GMP_MPF_BITS);
      mpf_set (v2, c.v2);
    }

    inline prec_t& operator= (const prec_t& c)
    {
      mpf_set (v2, c.v2);
      return *this;
    }

    inline prec_t& operator= (double v)
    {
      mpf_set_d (v2, v);
      return *this;
    }

    inline ~prec_t ()
    {
      mpf_clear (v2);
    }

    friend inline bool operator== (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) == 0; }
    friend inline bool operator!= (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) != 0; }
    friend inline bool operator< (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) < 0; }
    friend inline bool operator<= (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) <= 0; }
    friend inline bool operator> (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) > 0; }
    friend inline bool operator>= (const prec_t& l, const prec_t& r) { return mpf_cmp(l.v2, r.v2) >= 0; }

    inline prec_t& operator += (const prec_t& r) { mpf_add (v2, v2, r.v2); return *this; }
    inline prec_t& operator -= (const prec_t& r) { mpf_sub (v2, v2, r.v2); return *this; }
    inline prec_t& operator *= (const prec_t& r) { mpf_mul (v2, v2, r.v2); return *this; }
    inline prec_t& operator /= (const prec_t& r)
    {
      if ( mpf_cmp_d (r.v2, 0) != 0 )
        {
          mpf_div (v2, v2, r.v2);
        }
      else
        {
          mpf_set_d (v2, 0);
        }
      return *this;
    }

    inline prec_t operator- () const
    {
      prec_t ret;
      mpf_neg (ret.v2, v2);
      return ret;
    }

    inline friend prec_t operator+ (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      mpf_add (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator- (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      mpf_sub (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator* (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      mpf_mul (ret.v2, l.v2, r.v2);
      return ret;
    }

    inline friend prec_t operator/ (const prec_t& l, const prec_t& r)
    {
      prec_t ret;
      if ( mpf_cmp_d (r.v2, 0) != 0 )
        {
          mpf_div (ret.v2, l.v2, r.v2);
        }
      else
        {
          mpf_set_d (ret.v2, 0);
        }
      return ret;
    }

    inline prec_t sqrt () const
    {
      prec_t ret;
      mpf_sqrt (ret.v2, v2);
      return ret;
    }

    inline prec_t abs () const
    {
      prec_t ret;
      mpf_abs (ret.v2, v2);
      return ret;
    }

    inline double d() const { return mpf_get_d(v2); }

  private:
    mpf_t v2;
  };
  inline prec_t prec_t_sqrt (const prec_t& v) { return v.sqrt(); }
  inline double prec2double (const prec_t& p) { return p.d(); }
  inline size_t prec2size_t (const prec_t& p) { return (size_t)p.d(); }
  inline double cos (const prec_t& v) { return cos(v.d()); }
  inline double acos (const prec_t& v) { return acos(v.d()); }
  inline double sin (const prec_t& v) { return sin(v.d()); }
  inline double asin (const prec_t& v) { return asin(v.d()); }
  inline double tan (const prec_t& v) { return tan(v.d()); }
  inline double atan (const prec_t& v) { return atan(v.d()); }
  inline double floor (const prec_t& v) { return floor(v.d()); }
  inline double ceil (const prec_t& v) { return ceil(v.d()); }
  inline double round (const prec_t& v) { return round(v.d()); }
  inline double lround (const prec_t& v) { return lround(v.d()); }
  inline double prec_t_pow (const prec_t& v, const prec_t& f) { return pow(v.d(), f.d()); }
  inline double prec_t_log (const prec_t& v) { return log(v.d()); }
  inline prec_t fabs (const prec_t& v) { return v.abs(); }
#else
  typedef double prec_t;
  #define prec_t_sqrt sqrt
  #define prec_t_log log
  #define prec_t_pow pow
  #define prec2double(v) ((double)(v))
  #define prec2size_t(v) ((size_t)(v))
  #define prec_checkerror(a,b,c)
  #define prec_clearerror(a)
  #define prec_correcterror(a)
  #define PREC_DOUBLE
#endif

#if defined(_WIN32) || defined(_WINDOWS)
#define ATTRIBUTE_PACKED
#define ATTRIBUTE_ALIGNED16
// oh my god, why do I have to do this?
__forceinline double prec_t_round (double v)
{
    int result;
    __asm
    {
        fld    v
        fistp  result
    }
    return result;
}

__forceinline long int prec_t_lrand ()
{
  return rand();
}

__forceinline prec_t prec_t_drand ()
{
  return (prec_t)rand()/(prec_t)RAND_MAX;
}
#else
#define ATTRIBUTE_PACKED __attribute__ ((packed))
#define ATTRIBUTE_ALIGNED16 __attribute__ ((aligned (16)))
#define prec_t_round round
#define prec_t_lrand lrand48
#define prec_t_drand drand48
#endif

template <class T>
T prec_t_trunc (const T& v)
{
  if ( v < 0 )
      return -floor(-v);
  else
      return floor(v);
}

#define PREC_PI    3.141592653589793238462643383279502884197169399375105820974944592304
#define PREC_SQRT2 1.414213562373095048801688724209698078569671875376948073176679737990

/** Generic vector for 3D space. */
struct vect3d
{
  prec_t x, y, z;

  inline vect3d() { }
  inline vect3d(prec_t a, prec_t b, prec_t c) : x(a), y(b), z(c) { }

  inline void clear () { x = 0.0; y = 0.0; z = 0.0; }
  inline void add (const vect3d& a) { x += a.x; y += a.y; z += a.z; }
  inline void normalize () { *this /= length(); }
  /** @return Normalized vector. */
  inline vect3d normal () const { vect3d ret = *this; ret.normalize(); return ret; }
  /** @return \f$x^2+y^2+z^2\f$ */
  inline prec_t magnitude () const { return *this * *this; }
  /** @return \f$\sqrt{x^2+y^2+z^2}\f$ */
  inline prec_t length () const { return prec_t_sqrt (magnitude()); }

  inline void operator*= (prec_t f) { x *= f; y *= f; z *= f; }
  inline void operator+= (const vect3d& v) { add(v); }
  inline void operator-= (const vect3d& v) { add(-v); }
  inline void operator/= (prec_t f) { x /= f; y /= f; z /= f; }
  // inline void operator/= (prec_t f) { prec_t m = 1 / f; x *= m; y *= m; z *= m; }
  inline void operator= (prec_t f) { x = f; y = f; z = f; }

  friend inline vect3d operator* (prec_t f, const vect3d& v)
  { return vect3d(v.x * f, v.y * f, v.z * f); }

  friend inline vect3d operator* (const vect3d& v, prec_t f)
  { return vect3d(v.x * f, v.y * f, v.z * f); }

  friend inline vect3d operator/ (const vect3d& v, prec_t f)
  { return vect3d(v.x / f, v.y / f, v.z / f); }

  friend inline vect3d operator+ (const vect3d& l, const vect3d& r)
  { return vect3d(l.x + r.x, l.y + r.y, l.z + r.z); }

  friend inline vect3d operator- (const vect3d& l, const vect3d& r)
  { return vect3d(l.x - r.x, l.y - r.y, l.z - r.z); }

  /** Cross product. */
  friend inline vect3d operator% (const vect3d& l, const vect3d& r)
  { return vect3d(l.y * r.z - l.z * r.y, l.z * r.x - l.x * r.z, l.x * r.y - l.y * r.x); }

  /** Dot product. */
  friend inline prec_t operator* (const vect3d& l, const vect3d& r)
  { return l.x * r.x + l.y * r.y + l.z * r.z; }

  friend inline bool operator== (const vect3d& l, const vect3d& r)
  { return l.x == r.x && l.y == r.y && l.z == r.z; }

  friend inline bool operator!= (const vect3d& l, const vect3d& r)
  { return !(l == r); }

  friend inline bool operator< (const vect3d& l, const vect3d& r)
  {
    if ( l.x < r.x )
        return true;
    else if ( l.x == r.x && l.y < r.y )
        return true;
    else
        return l.x == r.x && l.y == r.y && l.z < r.z;
  }

  inline vect3d operator- () const { return -1.0 * *this; }

  inline vect3d recip () const { return vect3d (1.0 / x, 1.0 / y, 1.0 / z); }

  /** Find a vector that is at right angle to this. */
  void across (vect3d& out);
  
  /** Print vector value into stdout. */
  void print(const char *name = 0) const;
  /** Print vector into buffer. */
  char * sprint(char *buf) const;
};

/** 3D transformation matrix. No support for affine transformations. */
struct transf3d
{
  prec_t a, b, c;
  prec_t d, e, f;
  prec_t g, h, i;

  transf3d () { }

  transf3d (const vect3d& r1, const vect3d& r2, const vect3d& r3)
  {
    a = r1.x; b = r1.y; c = r1.z;
    d = r2.x; e = r2.y; f = r2.z;
    g = r3.x; h = r3.y; i = r3.z;
  }

  transf3d (prec_t f1, prec_t f2, prec_t f3, prec_t f4, prec_t f5, prec_t f6, prec_t f7, prec_t f8, prec_t f9)
  {
    a = f1; b = f2; c = f3;
    d = f4; e = f5; f = f6;
    g = f7; h = f8; i = f9;
  }

  /** Print transformation matrix into the stdout. */
  void print(const char *name = 0) const;

  /** Create an identity transformation. */
  void identity ()
  {
    a = e = i = 1;
    b = c = d = f = g = h = 0;
  }

  /** Create a reflection transformation based on a plane vector. */
  void reflect (const vect3d& n)
  {
    a = 1 - 2 * n.x * n.x;  b = -2 * n.x * n.y;     c = -2 * n.x * n.z;
    d = -2 * n.x * n.y;     e = 1 - 2 * n.y * n.y;  f = -2 * n.y * n.z;
    g = -2 * n.x * n.z;     h = -2 * n.y * n.z;     i = 1 - 2 * n.z * n.z;
  }

  /** Transform the given vector. */
  vect3d transform(const vect3d& v) const
  {
    return vect3d (a * v.x + b * v.y + c * v.z,
        d * v.x + e * v.y + f * v.z,
        g * v.x + h * v.y + i * v.z);
  }

  /** Calculate the determinant of this matrix. */
  prec_t det() const
  {
    return a * e * i + c * d * h + b * f * g -
        a * f * h - b * d * i - c * e * g;
  }

  /** Transpose of the matrix. */
  transf3d transpose () const
  {
    return transf3d (
        a, d, g,
        b, e, h,
        c, f, i);
  }

  /** Adjugate of the matrix. */
  transf3d adj () const
  {
    return transf3d (
        e*i-f*h,    -(b*i-c*h),  b*f-c*e,
        -(d*i-f*g), a*i-c*g,     -(a*f-c*d),
        d*h-e*g,    -(a*h-b*g),  a*e-b*d
    );
  }

  transf3d inverse () const
  {
    return adj() / det();
  }

  friend transf3d operator/ (const transf3d& t, prec_t f)
  {
    return transf3d (
        t.a / f, t.b / f, t.c / f,
        t.d / f, t.e / f, t.f / f,
        t.g / f, t.h / f, t.i / f);
  }

  friend transf3d operator* (const transf3d& t, prec_t f)
  {
    return transf3d (
        t.a*f, t.b*f, t.c*f,
        t.d*f, t.e*f, t.f*f,
        t.g*f, t.h*f, t.i*f);
  }

  friend transf3d operator* (prec_t f, const transf3d& t)
  {
    return transf3d (
        t.a*f, t.b*f, t.c*f,
        t.d*f, t.e*f, t.f*f,
        t.g*f, t.h*f, t.i*f);
  }

  friend vect3d operator* (const transf3d& t, const vect3d& v)
  {
    return t.transform(v);
  }

  friend transf3d operator * (const transf3d& l, const transf3d r)
  {
    transf3d ret;
    ret.a = l.a * r.a + l.b * r.d + l.c * r.g;
    ret.b = l.a * r.b + l.b * r.e + l.c * r.h;
    ret.c = l.a * r.c + l.b * r.f + l.c * r.i;
    ret.d = l.d * r.a + l.e * r.d + l.f * r.g;
    ret.e = l.d * r.b + l.e * r.e + l.f * r.h;
    ret.f = l.d * r.c + l.e * r.f + l.f * r.i;
    ret.g = l.g * r.a + l.h * r.d + l.i * r.g;
    ret.h = l.g * r.b + l.h * r.e + l.i * r.h;
    ret.i = l.g * r.c + l.h * r.f + l.i * r.i;
    return ret;
  }
};

struct transf3d_pair
{
  transf3d f;
  transf3d b;
};

/** Calculate complete elliptic integral of the first kind. */
extern prec_t ellipse1(prec_t kk);
/** Calculate complete elliptic integral of the second kind. */
extern prec_t ellipse2(prec_t kk);
/** Calculate complete elliptic integral of the second kind.
  * @param e1v Value of complete elliptic integral of the first kind for the same kk.
  */
extern prec_t ellipse2(prec_t kk, prec_t e1v);

extern prec_t prec_t_sub (prec_t a, prec_t b);

#endif // !__math3d_hpp__

